import json
import numpy as np
import pandas as pd
import scipy.optimize
import scipy.stats as st
import scipy.special
import numba
import tqdm
import biocircuits
import warnings
import iqplot
import bokeh.io
bokeh.io.output_notebook()
import panel as pn
pn.extension()
warnings.filterwarnings('ignore')
def style(p, autohide=False):
p.title.text_font="Helvetica"
p.title.text_font_size="16px"
p.title.align="center"
p.xaxis.axis_label_text_font="Helvetica"
p.yaxis.axis_label_text_font="Helvetica"
p.xaxis.axis_label_text_font_size="13px"
p.yaxis.axis_label_text_font_size="13px"
p.background_fill_alpha = 0
if autohide: p.toolbar.autohide=True
return p

$\hspace{15em}\mathrm{Nondimensionalizing... }\\[0.5em]$ \begin{align} \cfrac{c_d}{t_d}\cfrac{\mathrm{d}\tilde{c}}{\mathrm{d}t} &= \beta \cfrac{n_d\tilde{n}/k}{1+n_d\tilde{n}/k} c_d \tilde{c} - \cfrac{Q}{V}c_d \tilde{c} \\[0.5em] \cfrac{n_d}{t_d}\cfrac{\mathrm{d}\tilde{n}}{\mathrm{d}t} &= \cfrac{Qn_0}{V} - \cfrac{Qn_d\tilde{n}}{V} - \nu \beta \cfrac{n_d\tilde{n}/k}{1+n_d\tilde{n}/k} c_d \tilde{c} \\[2.5em] \end{align}
$\hspace{15em}\mathrm{Setting... }$ \begin{align} \boxed{n_d = k, \hspace{0.5em} t_d = \cfrac{V}{Q}} \\[1.0em] \end{align}
\begin{align} \cfrac{\mathrm{d}\tilde{c}}{\mathrm{d}t} &= \beta \cfrac{V}{Q} \cfrac{\tilde{n}}{1+\tilde{n}} \tilde{c} - \tilde{c} \\[0.5em] \cfrac{\mathrm{d}\tilde{n}}{\mathrm{d}t} &= \cfrac{n_o}{k} - \tilde{n} - \nu \beta \cfrac{V}{Q} \cfrac{c_d}{k} \cfrac{\tilde{n}}{1+\tilde{n}}\tilde{c} \\[0.5em] \end{align}$\hspace{15em}\mathrm{Setting... }$ \begin{align} \boxed{ \tilde{\beta} = \beta \cfrac{V}{Q}, \hspace{0.8em} \tilde{n_0} = \cfrac{n_0}{k}, \hspace{0.8em} c_d = k } \end{align}
\begin{align} \cfrac{\mathrm{d}\tilde{c}}{\mathrm{d}t} &= \tilde{\beta} \cfrac{\tilde{n}}{1+\tilde{n}} \tilde{c} - \tilde{c} \\[0.5em] \cfrac{\mathrm{d}\tilde{n}}{\mathrm{d}t} &= \tilde{n_0} - \tilde{n} - \nu \tilde{\beta} \cfrac{\tilde{n}}{1+\tilde{n}}\tilde{c} \end{align}$\hspace{15em}\mathrm{Rewriting... }$ \begin{align} \boxed{ \cfrac{\mathrm{d}c}{\mathrm{d}t} = \beta \cfrac{n}{1+n} c - c \\[0.5em] \cfrac{\mathrm{d}n}{\mathrm{d}t} = n_0 - n - \nu \beta \cfrac{n}{1+n}c } \end{align}
$\hspace{10em}\text{We find the nullclines and fixed points by setting the derivatives to zero:}$ \begin{align} \mathrm{nullcline}_1. \hspace{0.5em} &c = 0, && n = n_0 \\[0.5em] \mathrm{nullcline}_2. \hspace{0.5em} &c = \cfrac{(n_0 - n)(1+n)}{\nu \beta n}, && n = \cfrac{1}{\beta-1} \\[0.5em] \mathrm{fp}_1. \hspace{0.5em} &c_{st} = 0, &&n_{st} = n_0 \\[0.5em] \mathrm{fp}_2. \hspace{0.5em} &c_{st} = \cfrac{1}{\nu}\left(n_0 - \cfrac{1}{\beta-1} \right), && n_{st} = \cfrac{1}{\beta - 1} \\[0.5em] \end{align}
$\hspace{16em}\text{Note that for $\mathrm{fp}_2$ to be valid, $\boxed{\beta > 1, \hspace{0.5em} n_0 \geq \cfrac{1}{\beta-1}}$} $
$\hspace{10em}\text{We can then write out the linear stability matrix:}$
\begin{align} A = \begin{pmatrix} \beta \cfrac{n}{1+n} - 1 & \beta c \cfrac{1}{(1+n)^2} \\ -\nu \beta \cfrac{n}{1+n} & -1 -\nu \beta c \cfrac{1}{(1+n)^2} \end{pmatrix} \end{align}$\hspace{10em}\text{Using Mathematica, the eigenvalues are:}$
\begin{align} \lambda_i = -1, \hspace{2em} \lambda_{ii} = \cfrac{(\beta-1)n^2 + (\beta-2)n - (\beta c \nu + 1) }{(1+n)^2} \end{align}At the fixed point $\mathrm{fp}_1$ \begin{align} \lambda_{ii} &= \cfrac{n_0 (\beta-1) - 1}{1+n_0} < 0 \\[0.5em] &\text{when either 1. } \beta < 1,\text{or 2. } n_0 < \cfrac{1}{\beta-1} \hspace{0.1em}\mathrm{and}\hspace{0.1em} \beta > 1. \end{align}
At the fixed point $\mathrm{fp}_w$ \begin{align} \lambda_{ii} &= 1-(\beta-1)n_0 < 0 \\[0.5em] &\mathrm{if}\hspace{0.6em} n_0 > \cfrac{1}{\beta - 1}\hspace{12em} \end{align}
Note that for $\mathrm{fp}_2$ we take $\beta$ to be greater than 1 since its existence is predicated on this condition.
We can break up our analysis into two cases:
Production is not fast enough relative to how quickly $c$ is being flushed out. Only one fixed point exists where the bacterial concentration decays to zero, and the nutrient concentration steadies to $n_0$. This fixed point has two negative eigenvalues (look at numerator of $\lambda_{ii}$), and is thus stable.
When production is fast enough, we also need the influx nutrient concentration to be sufficiently high, $n_0 \geq 1/(\beta-1)$ to reach a nonzero bacterial concentration. When both conditions are satisfied, fixed point 2 (nonzero $c$) is stable, and fixed point 1 (zero $c$) is unstable (condition 2 is not satisfied since the inequality is flipped).
When $n_0$ is not sufficiently high, fixed point 1 (zero $c$) is again stable (condition 2 is satisfied once more).
First some colors.
color_C = "#696DAA"
color_N = "#C6A1B7"
color_phase = "#BCBCBC"
color_C_null = color_C
color_N_null = color_N
Define some functions next.
def dc_dt(c, n, β, n0, ν):
return β*n/(1+n)*c - c
def dn_dt(c, n, β, n0, ν):
return n0 - n - ν*β*n/(1+n)*c
def derivs(x, t, β, n0, ν):
c, n = x
dc_dt = β*n/(1+n)*c - c
dn_dt = n0 - n - ν*β*n/(1+n)*c
return np.array([dc_dt, dn_dt])
def get_fps(β, n0, ν):
'''returns (c_st, n_st)'''
fp1 = np.array([0, n0])
fp2 = np.array([1/ν*(n0 - 1/(β-1)), 1/(β-1)])
return fp1, fp2
def test_stability(args, fp):
β, n0, ν = args
c, n = fp
λ_numerator = (β - 1)*n**2 + (β-2)*n - (β*c*ν + 1)
stab = "stable" if λ_numerator < 0 else "unstable"
return stab
def draw_fp(x, y, p, fp_type="stable", in_radius=0.5, out_radius=0.2, n_wedge=12):
if fp_type == "stable":
p.circle(x, y, size=16, color="black", line_color='white', line_width=0.8)
if fp_type == "unstable":
p.circle(x, y, size=12, color="white", line_color="black", line_width=4)
# pizazz wedges
spin = 0
start_angles = np.linspace(0, 2*np.pi, n_wedge)
end_angles = start_angles + 0.1
in_radius = in_radius
out_radius = out_radius
for start_angle, end_angle in zip(start_angles, end_angles):
p.annular_wedge(
x=x,
y=y,
fill_color="black",
line_color="black",
inner_radius=in_radius,
outer_radius=in_radius+out_radius,
start_angle=start_angle+spin,
end_angle=end_angle+spin,
)
return p
def draw_phase(c_range, n_range, args, p, color="#B0B0B0", density=2):
# .... PHASE PORTRAIT ....
p = biocircuits.phase_portrait(
dc_dt,
dn_dt,
c_range,
n_range,
args,
args,
color=color_phase,
p=p,
density=2
)
return p
def draw_nullclines(c_range, n_range, args, p, color_C_null, color_N_null, lw=3, _shape=200, draw_zero=False):
β, n0, ν = args
c_zero = np.zeros(_shape)
n_n0 = np.ones(_shape) * n0
p.line(c_zero, (n_range[0]-2, n_range[1]+2), line_width=lw, color=color_C_null)
p.line((c_range[0]-2, c_range[1]+2), n_n0, line_width=lw, color=color_C_null)
c_space = np.linspace(c_range[0]-2, c_range[1]+2, _shape)
n_space = np.linspace(n_range[0]-2, n_range[1]+2, _shape)
c_n_null = (n0-n_space)*(1+n_space) / (ν*β) / n_space
n_n_null = np.ones(_shape) * (1/(β-1))
p.line(c_space, n_n_null, line_width=lw, color=color_N_null)
p.line(c_n_null, n_space, line_width=lw, color=color_N_null)
if draw_zero:
n_zero = np.zeros(_shape)
p.line(c_space, n_zero, line_width=lw*1.5, color="#d6543a", line_dash="solid")
p.line(c_zero, n_space, line_width=lw*1.5, color="#d6543a", line_dash="solid")
return p
def draw_legend(color_C_null, color_N_null, p, lw=3):
legend = bokeh.models.Legend(items=[
("stable fp", [p.circle(color="black", size=15)]),
("unstable fp", [p.circle(color="white", size=15, line_color="black", line_width=2)]),
("dc/dt nullclines", [p.line(color=color_C_null, line_width=lw)]),
("dn/dt nullcline", [p.line(color=color_N_null, line_width=lw)]),
], location="center")
p.add_layout(legend, 'right')
return p
def draw_everything(c_range, n_range, args,
color_C_null, color_N_null, color_phase,
lw=3.5,
density=2,
in_radius=0.15,
out_radius=0.05,
n_wedge=13,
title="Nullclines on the 𝑐−𝑛 plane",
draw_zero=False,
show_legend=True,
):
fp1, fp2 = get_fps(*args)
width = 600 if show_legend else 460
p = bokeh.plotting.figure(
height=400,
width=width,
title=title,
x_axis_label="c",
y_axis_label="n",
x_range=c_range,
y_range=n_range,
)
# .... PHASE PORTRAIT ....
p = draw_phase(c_range, n_range, args, p, color=color_phase, density=2)
# .... NULLCLINES ....
p = draw_nullclines(c_range, n_range, args, p, color_C_null, color_N_null, lw=lw, draw_zero=draw_zero)
# .... FIXED POINTS ....
stab1 = test_stability(args, fp1)
stab2 = test_stability(args, fp2)
p = draw_fp(fp1[0], fp1[1], p, fp_type=stab1, in_radius=in_radius, out_radius=out_radius, n_wedge=n_wedge)
p = draw_fp(fp2[0], fp2[1], p, fp_type=stab2, in_radius=in_radius, out_radius=out_radius, n_wedge=n_wedge)
# .... LEGEND ....
if show_legend:
p = draw_legend(color_C_null, color_N_null, p, lw=lw)
p = style(p, autohide=True)
return p
Let's call all our functions! I am setting $\beta > 1$ and $n_0 > 1/(\beta-1)$, and we see exactly what we expect, with the entire positive domains swirling into our steady state.
β, n0, ν = 1.2, 10, 5
args = (β, n0, ν)
c_range = (-0.5, 5)
n_range = (-0.5, 12)
params_c, params_m = args, args
fp1, fp2 = get_fps(*args)
lw = 3.5 # nullclines & legend
density = 2 # phase
in_radius = 0.16 # unstable fp
out_radius = 0.07 # unstable fp
n_wedge=15
p = draw_everything(
c_range, n_range, args, color_C_null, color_N_null, color_phase,
lw=lw, density=density, in_radius=in_radius, out_radius=out_radius, n_wedge=n_wedge,
title=f"Nullclines β:{np.round(β, 2)}"
)
bokeh.io.show(p)
Let's look at the three cases outlined above. On the left, we have the case when both conditions are satisfied. In the middle, n0 is not sufficiently large. On the right, $\beta$ is less than 1. In both cases we venture into the negatives.
c_range = (-1.5, 6)
n_range = (-6, 12)
in_radius = 0.2
out_radius = 0.08
plots = [
draw_everything(
c_range, n_range, (β, n0, ν), color_C_null, color_N_null, color_phase,
lw=lw, density=density, in_radius=in_radius, out_radius=out_radius, n_wedge=n_wedge,
title=f"Nullclines β:{np.round(β, 2)}, n0: {n0}", draw_zero=True, show_legend=show_legend
) for β, n0, show_legend in zip([1.2, 1.2, 0.8], [10, 2, 10], [False, False, True])
]
bokeh.io.show(bokeh.layouts.layout([plots]))
Note the subtle changes in the stabilities. I plot things in the negatives to get a better idea of what's happening with the math. Big red lines are drawn on the zeroes.
Let's plot some trajectories. Note that my initial conditions are thrown in randomly, but the random distribution is over a log scale to get better resolution at lower concentrations.
β, n0, ν = 3, 5, 2
args = (β, n0, ν)
seed = np.random.seed(12345789)
N_SEED, SCALE = 30, 2
Cos = np.random.random(N_SEED) * SCALE - 1.3
Nos = np.random.random(N_SEED) * SCALE - 1.3
Cos, Nos = np.power(10.0, Cos), np.power(10.0, Nos)
β_slider = pn.widgets.FloatSlider(name="β", start=0.15, end=5, value=1.85, step=0.1, width=170)
n0_slider = pn.widgets.FloatSlider(name="n0", start=2, end=12, value=9, step=1, width=170)
ν_slider = pn.widgets.FloatSlider(name="ν", start=1.2, end=5, value=3.3, step=0.1, width=170)
t_max_slider = pn.widgets.IntSlider(name="t max", start=5, end=30, value=12, step=5, width=300)
@pn.depends(β_slider.param.value, n0_slider.param.value, ν_slider.param.value, t_max_slider.param.value)
def plotter_traj(β, n0, ν, t_max):
args = (β, n0, ν)
t = np.linspace(0, t_max, 100)
q = bokeh.plotting.figure(
height=450, width=550,
title='Trajectories',
x_axis_label='dimensionless t',
y_axis_label='[ ]',
y_range=(-0.2, 9)
)
for C, N in zip(Cos, Nos):
CNo = np.array([C, N])
_CN = scipy.integrate.odeint(derivs, CNo, t, args=args)
C, N = _CN.T
q.line(t, N, line_color=color_N, line_width=2.7, line_alpha=0.5, legend_label="n")
q.line(t, C, line_color=color_C, line_width=2.7, line_alpha=0.5, legend_label="c")
if (n0 >= 1/(β-1)) and (β > 1):
c_st = c_st = 1/ν * (n0 - 1/(β-1))
n_st = 1/(β-1)
else:
c_st = 0
n_st = n0
q.line((0, t.max()), 0, line_color="black", line_width=3.0, line_dash="dashed", legend_label="0")
q.line((0, t.max()), c_st, line_color="#696D9A", line_width=5, line_dash="dotdash", legend_label="c_st")
q.line((0, t.max()), n_st, line_color="#B691C7", line_width=5, line_dash="dotdash", legend_label="n_st")
q.legend.location = "top_right"
return style(q, autohide=True)
lay_widgets = pn.Column(β_slider, n0_slider, ν_slider, align="center")
lay_time = pn.Row(t_max_slider, align="center")
dashboard = pn.Row(pn.Column(plotter_traj, lay_time), lay_widgets)
dashboard
Here's a quick snapshot of the different behaviors.

Toggling $\beta$ and $n_0$, we get precisely the trend we expect, where $c$ falls in the regime of low $\beta < 1$, and rises to a steady state when both n0 and β are large. Increasing nu decreases our steady state $c$, and increasing n0 increases our steady state $c$. Since nu is harder to tune in real life, what we have to do to experimentally is:
The above inequalities yield the following. \begin{align} Q & < \beta V \\[0.5em] Q & < \beta V \cfrac{1}{1 + k/n_0} \end{align} Note that the second condition yields the stricter inequality. Thus the chosen inlet flow rate should be low enough to give the system enough time to build up sufficient quantities of $c$ before it gets flushed out. $$\boxed{Q < \beta V \cfrac{1}{1 + k/n_0}}$$
Since $\beta$, $k$, and $\nu$ are intrinsic to the nutrients and bacteria, we can tune $n_0$ to give the desired steady-state $c$.
Note: we would have to solve Q and n_0 simultaneously. The general idea in practice is to first try a pair of ($n_0$, Q) with small Q/V relative to $\beta$, and $n_0$ close to $k$ to then yield a desirable $\tilde{c_{st}}$, then adjust Q according to the final inequality. For some experimental values of $V, k, \beta, \nu$, we can also computationally find the perfect pair.
import perl
palette = perl.palette
k = 1
β = 1
V = 1
ν = 1
# *************** NOTE a: no, b: Q ***************
na, nb = 150, 150
a = np.linspace(0.1, 10.0, na)
b = np.logspace(-3, -0.1, nb)
aa, bb = np.meshgrid(a, b, indexing='xy')
C_grid = np.empty((na, nb))
for i in range(na):
for j in range(nb):
_a, _b = aa[i, j], bb[i, j]
C_grid[i, j] = 1/ν * (_a - k/(β*V/_b-1))
_df_alphas, _df_betas, _df_metric = [], [], []
for i in range(na):
for j in range(nb):
alpha, beta = aa[i, j], bb[i, j]
metric = C_grid[i, j]
_df_alphas.append(alpha)
_df_betas.append(beta)
_df_metric.append(metric)
_df = pd.DataFrame({'no': _df_alphas, 'Q': _df_betas, 'c_st':_df_metric})
q = bokeh.plotting.figure(
title="steady [c] in 𝑄-no plane",
width=450, height=400,
y_axis_type="log",
x_axis_label="no",
y_axis_label="𝑄"
)
colormap = bokeh.models.LinearColorMapper(palette=palette)
colorbar = bokeh.models.ColorBar(color_mapper=colormap, location=(0,0))
q.circle(source=_df, x='no',y='Q', color={'field':'c_st', 'transform':colormap}, size=4)
no_space = np.linspace(a[0], a[-1], 1000)
q.line(no_space, β * V / (1+k/no_space), line_width=10, line_color="grey")
q.line(no_space, β * V / (1+k/no_space), line_width=7, line_color="black")
q.add_layout(colorbar, 'right')
bokeh.io.show(style(q))
The color represents the [c]. The black line represents $Q = \beta V / (1+k/n_0)$. All points below this curve satisfy both relations. We see that above the curve, [c] dips into the negatives. We thus see that we can pass in the fixed values of all the other parameters, locate the desired $c$ on the color-bar, then pick an ($n_0$, $Q$) pair. $Q$ will partly control the timescales of the experiment (how long it takes to reach this steady state), so the fact we have a log-scale to choose from is very nice.
In this problem, I really enjoyed how elegant the results were even though it was a relatively simple system, because I felt like I understood why every parameter would be important in the context of an experiment. Chemostat, thank you for the chat! And reader, as always, thank you for reading!